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ABSTRACT 

We examine the structure of the post-shock region in supernova remnants 
(SNRs). The "shock transition zone" is set up by charge transfer and ionization 
events between atoms and ions, and has a width ~ 10^^ cm~^ tt-q^, where uq 
is the total pre-shock density (including both atoms and ions). For Balmer- 
dominated SNRs with shock velocity Vg > 1000 km s~^, the Rankine-Hugoniot 
conditions for ion velocity and temperature are obeyed instantly, leaving the 
full width at half-maximum (FWHM) of the broad Ha line versus Vg relation 
intact. However, the spatial variation in the post-shock densities is relevant to 
the problem of Lya resonant scattering in young, core-collapse SNRs. Both 
two- (pre-shock atoms and ions) and three-component (pre-shock atoms, broad 
neutrals and ions) models are considered. We compute the spatial emissivities 
of the broad (^f,) and narrow (^„) Ha lines; a calculation of these emissivities 
in SN 1006 is in general agreement with the computed ones of Raymond et 
al. (2007). The (dimensionless) spatial shift, 9shift, between the centroids of 
and C,n is unique for a given shock velocity and /ion, the pre-shock ion fraction. 
Measurements of Gshift can be used to constrain uq. 

Subject headings: atomic processes — hydrodynamics — shock waves — super- 
nova remnants 



1. INTRODUCTION 

In a purely hydrodynamical treatment, shock fronts are regarded as mathematical dis- 
continuities, across which the density, pressure and temperature of the fluid vary according to 
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the Rankine-Hugoniot jump conditions. A characteristic length scale for the thickness of the 
discontinuity appears when one takes into account the atomic structure of the gas (Zel'dovich 
&; Raizer 1966). The structure of coUisional shock fronts is relatively well understood. A 
dissipation mechanism is generally required: for example, thermal conduction and viscosity 
for weak and strong shocks, respectively. In collisionless shock fronts, the dissipation is due 
to turbulent electromagnetic fields rather than collisions. Determining the structure of col- 
lisionless shocks requires an understanding of how particles equilibrate their temperatures 
via plasma instabilities, which is currently a largely open question in astrophysics. Reviews 
can be found in McKee & Hollenbach (1980) and Draine & McKee (1993). 

Our study is motivated by a need to understand the structure of the post-shock region 
in Balmer-dominated supernova remnants (SNRs; Heng & McCray 2007, hereafter HM07, 
and references therein), which we term the "shock transition zone". These SNRs typically 
have shock velocities Vs ~ 1000 km s~^, ages much less than the radiative cooling times 
and produce a modest amount of ionizing radiation. Balmer-dominated SNRs have two- 
component spectra consisting of a broad and a narrow line (Chevalier, Kirshner & Raymond 
(1980, hereafter CKR80; see §6.11) . If the temperatures of the atoms and ions are known, the 
width of the broad line is uniquely related to the shock velocity, as first shown by CKR80 
The upstream, pure hydrogen gas consists of atoms and ions (T^ ~ 10000 K) with a typical 
pre-shock ion fraction of /ion ~ 0.5. The atoms are converted into ions via charge transfer 
(with protons) and impact ionization (with both electrons and protons) in the transition 
zone (Fig. [T]), which has a width on the order of the mean free path of interactions. 

If a noticeable variation of the ion velocity occurs within the transition zone, then for a 
given broad line width, the shock velocity would be under- estimated. This is because charge 
transfer events, which give rise to broad Ha emission, are favored at low relative velocities 
between the atoms and ions, which may occur within the transition zone. In §6.H we show 
this not to be the case. 

The second reason for investigating the structure of Balmer-dominated SNRs is that 
Lya resonant scattering occurs in the freely streaming debris of young, core-collapse SNRs 
that are still in the pre- Sedov- Taylor phase. Photons produced at the shock fronts undergo 
a spectral random walk and become increasingly redshifted as the debris is in Hubble-like 
flow. Lya line profiles will therefore be distorted with respect to non-resonant lines such 
as Ha (if the hydrogen atoms are mostly in the ground state). Photon production occurs 
within the shock transition zone, the width of which is usually greater than the mean free 
path of the Lya photons. Therefore, understanding the spatial structure of the transition 
zone in Balmer-dominated SNRs is relevant to modeling the Lya lines in young remnants 
like SNR 1987A, where resonant scattering occurs in the hydrogen ejecta from the massive 
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progenitor. 

The third motivation for our study is to develop a new method for deriving the spatial 
emissivity profiles of both the broad and narrow Ha lines behind the shock front. Raymond 
et al. (2007, hereafter R07) have shown that such profiles can be used to infer no, the total 
pre-shock density (including both atoms and ions). 

In §21 we present our model and assumptions and state the relevant equations involved. 
We describe our solution methods in §1] and analyze our results in §3 In ^ we discuss the 
implications and limitations of our results and present opportunities for future work. 



2. MODEL & ASSUMPTIONS: TWO-COMPONENT MODEL 

For the problem to be tractable, we need to make a few simplifying assumptions. Firstly, 
the only significant sink term present for the atoms is for their conversion to ions, via charge 
transfer and impact ionization. In the frame of the shock, the atoms comprise a beam with 
a velocity of vu = Vg- The ions are shocked and isotropized in a distance comparable to 
the cyclotron gyroradius (/gyro ~ 10^° cm). We assume that the ion distribution function 
becomes approximately Maxwellian in a length scale that is much less than the width of 
the shock transition zone. Plasma instabilities are capable of partially equilibrating the 
temperatures between the electrons and protons (Te/Tp > 0.1 for Vg < 1000 km s~^; Cargill & 
Papadopoulos 1988; Rakowski, Ghavamian & Hughes 2003; Ghavamian, Laming & Rakowski 
2007). We parametrize the ratio of the temperatures at the shock front by 

P=^- (1) 

Hereafter, we drop the "p" subscript when referring to variables describing the ion fluid. The 
length scale for equilibration is determined by Coulomb collisions and can be much larger 
than the physical extent of the SNR for (3 = 0.1. It is generally larger than /gyro and the 
length scale for atomic interactions, /zone ~ 10^^ cm~^ Uq^. 



2.1. INTERACTION OF ATOMIC BEAM WITH MAXWELLIAN IONS 

In the shock transition zone, the atomic beam (number density nn and velocity fn) 
interacts with Maxwellian ions (density n, velocity v and temperature T). Charge transfers 
between the atomic beam and post-shock protons produce populations of atoms having 
velocity distributions intermediate between a beam and a Maxwellian, which we term "broad 
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neutrals"0. The resulting "skewed Maxwellian" distributions are described in HM07. These 
broad neutrals are produced in the same transition zone as the ions. 

We shall show how to account for the broad neutrals in §31 In this section, we assume 
low neutral fractions (/ion ^ 0.9) so that the broad neutrals will not affect the density of the 
ions significantly. In this case, we can approximate the variation of mass, momentum and 
enthalpy flux with distance by 

^('^h^^h) = -nunR, 



dx 

(P 

dx 



± {Uv + Pv + \pv^) = \mnvlnnnR, (2) 

where p = rripn. The internal energy density and pressure of the ion fluid are represented 
by U and P, respectively. The set of equations is supplemented by the following equations 
of state: 



P=(P + l) nkT, 

P ^ 3 
7-1 2 



U=^ = U(3 + l)nkT. (3) 



In the case of full equilibration {j3 = 1), we have P = 2nkT, where the factor of 2 accounts 
for contributions from both electrons and protons. The polytropic index is 7 = 5/3 for a 
monoatomic gas. The rate coefficient, R, characterizes the loss of atoms by conversion to 
ions via both charge transfer and ionization. Examples of R are shown in the left panel of 
Fig. [21 these rates assume a relative velocity of 3uh/4 between the atoms and the peak of 
the ion distribution. For /3 = 1, charge transfer dominates over ionization at vn < 4000 km 
s~^. Details on how to calculate R for a general relative velocity are given in Appendix 1X1 

Suppose we can define a typical value for the rate coefficient, R. Then a natural length 
scale for the problem emerges: 

L = ^. (4) 

UaR 

where is the pre-shock atomic density. Knowledge of L allows us to define the dimen- 
sionless distance, ( = x/L. Other dimensionless variables follow naturally: 77 = n/na, 



^In this paper, we refer to the pre-shock atoms in a beam simply as "atoms", while "broad neutrals" refer 
to the post-shock atoms found in broad distributions. Alternatively, the "atoms" could have been named 
"narrow neutrals" . We have avoided use of the terms "slow" and "fast neutrals" as these are dependent on 
one's frame of reference. 
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rjn = nfi/ua, u = f/fH, e = kT/rripV^, and TZ = R/R. The system of equations then 
becomes: 

^ [r]u^ + (/? + !) r^e] = r]nvn, 
^ [riu^ + 5 (/? + 1) rieu] = riuriTZ. (5) 



2.2. GENERALIZED CONSERVATION EQUATIONS 

The upstream, initial values for r/n, u and e are 1, = /ion/(l — /ion), 1 and 
eu = kTu/nipV^, respectively. Adding the first pair of equations in ^ and applying initial 
conditions, one obtains 

r]fi + r]u = l + fu, (6) 
which is a statement of the conservation of mass density flux, from which one can derive 

-T-r = (1 + /« - ^h) • (7) 

dC u 

Algebraic manipulation of the equations in ([5]) yields the cubic equation for u: 

{u - 1) [At]u^ -r]u^5iP+l) Uu] = 0, (8) 
With the solutions for rju, rj and u in hand, one can then solve for e using 

In general, since the equations in ([5]) have the same source/sink term, we can add/subtract 
any given pair and obtain algebraic equations equated to six conservation constants. Only 
three of the conservation constants are independent; the other three are a simple, linear com- 
bination of the first three. One can choose to tackle the system of four coupled differential 
equations or substitute up to three of them by the conservation equations. 



2.3. ENERGY LOSSES FROM THE ATOMIC BEAM 



The atomic beam suffers energy losses due to excitation and ionization, prior to engaging 
in charge transfer events with the ions. On average, an atom gets excited Ne ~ 1.5 times 
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before it gets taken out of the beam (HM07). Therefore, the final equation in ([5]) has to be 
modified: 

— [r]u^ + 5 {(3+1) r]eu] = r^^r^TZ (1 - 2eio,,) . (10) 

The quantity eioss = -E'ioss/'^^h'^^h characterizes the energy loss, 

Eio,, = 13.6 eV + Ne (10.2 eV) . (11) 
The cubic equation for u changes slightly: 

4r/n=^ - hr^u^ + n [77 (1 - 2eioss) - 5 (/3 + 1) /.ej + 5 (/3 + 1) /.e, + 2eioss = 0. (12) 

We emphasize that this is an approximate way to account for the energy losses. More 
generally, Ne has a velocity dependence. However, eioss ~ 10~^ at fH ~ 1000 km s~^, and a 
more complicated treatment (e.g.. Cox & Raymond 1985) is only important for shocks with 
fH ^ 200 km s~^. Excitation and subsequent emission of Lyo; photons extract increments 
of 10.2 eV from the atomic beam; energy losses due to Balmer and other Lyman lines are a 
~ 10% effect for ~ 1000 km s"^ shocks. 



3. RESTRICTED THREE-COMPONENT MODEL 



At low and intermediate values of /ion, the creation of the broad neutrals competes effec- 
tively with that of the ions. In this case we must generalize the hydrodynamical equations to 
a three-component model, consisting of the beam of atoms ("H"), the electrons and protons 
("i"), and the broad neutrals ("B"). Here we assume that the broad neutrals and ions share 
a common fluid velocity, v-q = v, which implies that they share a common temperature as 
well (Tb = T). This approximation has been employed in all of the papers before HM07, 
who showed that the steep velocity dependence of the charge transfer cross section creates 
significant differences between f b and v (and hence Tb and T) for shock velocities > 3000 
km s~^. The full three-component model requires vb 7^ v, which is beyond the scope of our 
present study (see §6.2p . 

We consider three processes between the particles: charge transfer between the atoms 
and ions (with a rate coefficient of -Rch), ionization of atoms by ions {Rm), and ionization 
of broad neutrals by ions (RiB)] the rate coefficients are calculated using the formalism of 
HM07 and shown in the right panel of Fig. [2l We neglect interactions between the atomic 
beam and the broad neutrals. 

Under these assumptions, the equations in ([2]) generalize to: 

£ (nnvu) = -niin{Rc}i + Rm), 
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— (ubv) = riiinRcH - UbuRib, 
A (^ny^ = riunRai + UBuRiB, 
£ [{P + Pb) + \mpnv^ + imnnBt;^] = mnt^H'^H^ (^ch + Rm) , 
± [{U + Ub)v + {P + Pb) V + \mpnv^ + {mY^UBV^] = ^mnvl^nnn [Rcu + Rm) , (13) 

where P and U have the same definitions as before, while Pb = 't-b^Tb and Ub = 3Pb/2. 
Casting the equations in dimensionless form yields: 

^ = -r/H^(l + /iH)7^, 
^ (Vbu) = VnVT^ - VBVfiB^, 

I [{r]t + f37])e + r]tu^] = mv (1 + /m) 7^, 
I [5 (r^t + p7]) eu + r^^u^] = r^^V (1 + /ih) 7^, (14) 

where ?7h, ''7, and e retain their previous definitions, while rjB = TiB/ua, fij = Rij/Rcn, 
71 = Rcn/R and rjt = r] + rjB- The initial conditions are: ?7h(0) = 1, ''^(O) = /„, ?7b(0) = 0, 
u{0) = 1 and e(0) = e„. 

Our approach to solving for the hydrodynamical variables is similar to the one previously 
described in §2.2[ In fact, the cubic equation for u remains the same as equation ([8]), but 
with T] replaced by rjt: 

(u - 1) [Ar]tu^ -^^u-5{/3 + 1) /„e„] = 0. (15) 
Adding the first three equations in fll4p and applying initial conditions yields: 

+ = 1 + fu- (16) 
Defining ji = rju and /ib = ?7bw, we use equation (IT^ to eliminate ?7h and obtain: 

I = i (1 + /J /iH7^/i - ifmn^i^ - i {f,B - ifm) 7^/iB/i, 

Once rj and t]b are known, e is determined using 



^ = i (1 + A) 7^/i - i7^/i2 - i (1 + i/iB) 7^/iB/i. (i7) 



r^fU (1 - n) + (1 + /3) fueu 
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4. SOLUTION METHODS 

4.1. APPROXIMATE SOLUTIONS AND THEIR ASYMPTOTES FOR THE 

TWO-COMPONENT MODEL 

For t>H > 1000 km s~^, we can assume eioss ~ and constant IZ (see [|5]) behind the shock, 
and derive approximate solutions to the hydrodynamic variables and their corresponding 
asymptotes. Again assuming upstream initial conditions, we obtain from equation ([H]): 



u 



l± Wl + 



^0 (/3+l)/„e„ 

7] 



(19) 



We pick the positive root to ensure that m > 0. Consider the limiting case where ^ 1. 
We then have 

1 , 5(/? + l)Ae^_ ^20) 



M ^ - + 

4 



Using equation ([H]), one gets 



r/ ^ 4 [1 + /„ - r^H - 5 (/3 + 1) U^] . 



(21) 



Substituting this into the first equation in ^ yields the Bernoulli equation, r/n+aur/H = &u^H' 
where a„ = 47^ [1 + — 5(/3 + 1 )/„£„] and hu = 47?.. Solving for r/n yields: 



— + ( 1 - — ) exp (a„C) 



(22) 



The set of solutions, (m, t^h, e), has the following asymptotes (C ^ 1): 



u 


1 

4' 


7] 


-4, 


Vn 






3 




16{/3+l) ' 



(23) 

consistent with standard jump conditions for a strong shock. Another quantity of interest 
is the Mach number, M, of the post-shock ion fiow. Since P oc p^^^, we have = dP/dp = 
5(/3 + l)ef^/3. The Mach number then becomes 



(24) 
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4.2. NUMERICAL SOLUTIONS 

To obtain numerical solutions for the set of equations in ([5]), we employ a coordinate 
system where the beginning of the transition zone (Fig. [1]) is placed at ( = 0. We discuss 
our approach for eioss = 0; our solution method for eioss 7^ is conceptually similar. Firstly, 
we assume ?7h(0) = 1, and use equations (jH]), dHD, and Qj to solve for ^^(O), m(0), and e(0). 
To obtain a self-consistent set of solutions, we use an iterative approach. For the first 
iteration, the (dimensionless) ion velocity and temperature are taken to be u^^\() = u{0) 
and e*^^^(C) = e(0), respectively. Equation ([7]) is then solved using a standard Runge-Kutta 
method for ri^\C), which is substituted into equation ([6]) to determine ^^''^•'(C)- With ^^^^■'(C) 
known, the cubic equation in ([H]) is solved at each ( using a simple bisection method in the 
range < m < 1, for which there is only one physical root. (Note that if eioss <^ 1, we only 
need to solve a quadratic equation for u^'^\().) The updated value of u^'^\() is substituted 
into equation ([9]), yielding e^^^(C). This procedure is iterated i times, until the dependent 
variables VniC), v'^^iC), u'^'KC) and eW(C) converge. The convergence is monitored in two 
ways: 1. The final, fractional correction for each dependent variable must be less than a 
pre-determined tolerance, i.e., \u^'^\C) ~ '"'■*^^''(C)I/'^*"*'*(C) ^ ^toi; 2. The mass, momentum, 
and enthalpy fluxes must be constant and equal to their values at ( = 0. 

In practice, the values of the variables converge rapidly during the iteration, requiring 
i < 10. For shock velocities Vs ^ 1000 km s~^, the velocity difference between the atomic and 
ion populations is nearly constant, meaning R R {JZ ^ 1) throughout the shock transition 
zone. In this regime, eioss is negligible and the analytical solutions of §4.11 are excellent 
approximations to the two-component numerical calculations. At lower shock velocities 
(fH ^ 200 km s~^), variations in u{Q and non-neghgible values of eioss quantitatively change 
the solution, requiring the full numerical treatment. 

The numerical method in the three-component case is analogous to the iterative proce- 
dure described above. Assuming ?7h(0) = 1 and 77b(0) = 0, equations (fT5l) . (fT6l) . and (fTSjl are 
solved for ?7(0), m(0), and e(0). As before, we set m*^^^ = m(0), and now solve the equations in 
( TT7|) for n^^^ and /ig ■* using a standard Runge-Kutta algorithm, from which we calculate rj^^^ 
and ?7b^ using u'^^\ Equation f|T5l) is solved using the updated values for the (dimensionless) 
densities to give the improved estimate, u^'^\ The process is then iterated in the manner 
described above. 
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5. RESULTS 

Figure [3] displays the structure of the shock transition zone for vg = wh/IOOO km s~^ 
= 1, 5, 7 and 10. In each case, we assume /„ = 1 (/ion = 0.5) and T4 = T„/10^ K = 1. 
Knowledge of fg and T4 then determines e„ = 8.25 x 10^^ T^^/vl- The pre-shock atomic 
density is ria = n^lil + /„) ~ 0.1 cm"'^. For example, in the case of SN 1006, R07 find 
0.25 < no < 0.4 cm"^, implying 0.025 < < 0.04 cm"^. 

Hydrodynamical quantities vary from their pre- to post-shock values, at the beginning 
of the shock transition zone, according to the Rankine-Hugoniot jump conditions. Across the 
zone, the ions have no velocity structure, consistent with the assumption made by CKR80. 
Numerically, for a strong shock, u goes from 1 to 1/4 immediately and r] jumps to 4/^ to 
conserve ion mass flux. The latter is in accordance with the Rankine-Hugoniot density jump 
(Zel'dovich & Raizer 1966) of 

J _ (7 + 1) e^o + (7 - 1) e« 

'°"(7-l)Qo + (7 + l)e.' ^''^ 

where e^p is the value of e immediately after the shock (and not far downstream). For 
and 7 = 5/3, we recover the familiar density jump of 4 for a monoatomic gas. For 
strong shocks (fn ^ 1000 km s~^), energy losses from the atomic beam are negligible after 
the jump, and the downstream density eventually evolves to Jd ~ 4(1 + /«) ~ 8 for = 1. 
A factor of 4 comes from the jump condition, while an additional factor of 2 results from 
adding the atoms to the population of ions. We note that for weak shocks, a departure from 
the decrease by a factor of 4 in m occurs, consistent with equation (1251) . The departure from 
a jump of 4 in (and its subsequent evolution) follows naturally to conserve momentum; 
the asymptotic value of e dips to below 3/32 (/5 = 1) due to energy losses. It is worth noting 
that as the shock velocity increases, the distinction between 773 for /3 = 0.1 and 1 vanishes. 

With a telescope having sufficient angular resolution, it may be possible to measure the 
spatial emissivity profiles of the narrow and broad Ha lines. The narrow emissivity is given 
by = nn-aRfia,!!, where Rua,n is the rate coefficient for the narrow Ha line (cf. equation 
[20] in HM07). Broad emission has two distinct contributions: (1) from charge transfers of 
the atoms in the initial beam directly to excited states of the broad neutrals, with a rate 
coefficient given by Rna,ho] and (2) from excitations of and charge transfers between the 
broad neutrals to their excited states, with a net rate coefficient of -Rna.b*- The addition of 
-Rna.bo aiid -RHa,b, yiclds the broad-line rate coefficient (cf. equation [22] in HM07). Charge 
transfers involving atoms from the beam will naturally have the spatial profile of the narrow 
emissivity. 

We first derive the emissivity profile of the broad line, C,b, from using the two- 
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component model: upon its creation at xq, a broad neutral drifts for an average distance 
of Id = ld{x, Xq) until it gets destroyed by impact ionization. The drift velocity, Vd, is the 
velocity difference between the peak of the broad neutral distribution and the shock front, 
The broad Ha emissivity is then 



U^) = ^^^n + I Ux) P {x, x') dx\ (26) 

where P(x, a;o) = Pq{xq) exp [—{x — xo)/ld{x, xq)] is the "transfer function" . Details regard- 
ing ld{x, Xq) and P{x, Xq) are given in Appendix [Bl Examples of -Rho.ii, RRa,ho and Rua,h, 
are given in Fig. HI In the three-component model, the broad and narrow Ha emissivities 
are simply given by C,b = ?T,nH-RH«,bo + nnBR}ia,h, and ^„ = nnH-RHa,n, respectively. 

In the two-component model, the ratio of broad to narrow Ha emission is given by 

J ,j _ lo ^b{x) dx _ -RHa.bp + -RHa,b« 
Jo Cn[X) dx R}la,n 

such it is equal to the ratio of rate coefficients and the integrated line ratio is preserved. In 
the three-component model, this is not necessarily the 

-Rna.bo /o°° nnu dx + Rua,h, /o°° '^^b dx 

^bn — FBo J Th/in, (28) 

-KHa,n Jq '^^H dx 

due to the fact that Jq°° rann dx 7^ Jp°° nn-Q dx in general. 

The high shock velocity and low neutral fraction of the remnant of supernova (SN) 1006, 
located 2.1 kpc away, makes it an ideal case study for both the two- and three-component 
models. Following R07, we compute ^b and ^„ (Fig. [5]). We adopt the following parameters: 
wh = 2890 km s'^ (R07), Vd/vn = 0.34, /5 = 0.1 (Ghavamian et al. 2002, hereafter G02), 
/ion = fu/{l + fu) = 0.9 (G02), T4 = 1, = 0.025 cm'^ (no = 0.25 cm'^; R07), R = 
1.1 X 10-^ cm3 s-\ RcR = 5.9 x 10"^ cm^ s-\ Riu = 5.0 x 10"^ cm^ s'^ and Rib = 4.5 x 10"^ 
cm^ s^^. The broad Ha line has the following parameters: -Rno.bo = 1-5 x 10~^ cm^ s~^ and 
Rna,h, = 1-0 X 10"*^ cm^ s~^ (Case A conditions). The narrow Ha line has: Rna,n = 6.9 x 10~^ 
cm^ s~^ (Case A) and 1.7 x 10"^ cm^ s~^ (Case B). We emphasize that it is not our intention 
to model the Lyman line trapping, as done by Ghavamian et al. (2001) and G02. Rather, 
we wish to calculate the relative shifts between and and demonstrate that Lbn 7^ h/ In] 
these are not dependent on the opacity assumptions for the narrow Ha line. As a matter of 
illustration, we adopt Case B conditions for the narrow Ha line. 



Except for T4, values of the parameters without references were computed using the 
formalism of HM07. The drift velocity is f^/'i^H > 1/4 because at the velocity of SN 1006, 
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charge transfer is not efficient enougfi to create a Maxwellian population of broad neutrals 
centered at i'H/4. 

Our calculation shows that peaks at ~ (y/07 from the shock front, a factor of ~ 2 
smaller than the ~ 0'/14 value computed by R07. The smaller spatial scale of the current 
calculation is due to a numerical error in the model code used in R07; it implies that the 
lower values of no in the curved shock models of Fig. 5 of R07 will not produce too much 
Ha emission towards the inside of the remnant, and the density range 0.25 < Uq < 0.4 cm~^ 
obtained by R07 produces too steep a falloff. However, smaller pre-shock densities require 
larger radii of curvature to match the observed peak surface brightness, and that produces 
too gradual a falloff towards the outside of the remnant. Overall, the revised models of R07 
are compatible with densities in the range 0.15 < no < 0.3 cm~^, but with small-scale ripples 
in the shock front that broaden the filament by about 0'/5. The values of our emissivities are 
comparable to those of R07; minor discrepancies in the emissivities may be partially due to 
our use of cascade matrices to compute the Ha rate coefficients, following HM07. Excitation 
from the ground state tends to populate the lower / levels, especially the p ones. Hence, 
our calculations over-estimate the cascade contribution, while R07 ignores it, and the true 
emissivities are probably bracketed by these two approaches. 

For SN 1006, we compute h/In = 0.68 and Lhn = 0.77 (Case B), with the latter being 
about 13% higher than the former. Despite this difference, it is worthwhile to note that the 
shift between the peaks of ,^5 and is about the same in both the two- and three-component 
calculations. Furthermore, Smith ct al. (1991) and G02 measure the broad-to-narrow Ha 
line ratio to be about 0.73 and 0.84, respectively, and thus our prediction is within the range 
of uncertainty. 

6. DISCUSSION 

6.1. SNR 1987A & BALMER-DOMINATED SUPERNOVA REMNANTS 

Balmer-dominated SNRs are named for the dominance of their hydrogen fines (over 
forbidden ones), as first described by Chevafier & Raymond (1978) and CKR80. They are 
characterized by their two-component spectra, which consists of a narrow (~ 10 km s~^) 
superimposed upon a broad (~ 1000 km s~^) line. Broad line emission is produced when 
the atoms engage in charge transfer reactions with the post-shock ions. Narrow lines are the 
result of direct excitation of the pre-shock atoms. For a given value of /?, the full width at 
half-maximum (FWHM) of the broad component is uniquely related to the shock velocity, 
Vs = vh, a relation which provides a way of measuring the distances to SNRs (Kirshner, 
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Winkler & Chevalier 1987). 

As mentioned in §21 HM07 showed that the charge transfer cross section is sensitive to 
the shock velocity. For very fast shocks, the bulk velocity of the broad neutrals exceeds that 
of the protons, resulting in lower values of the FWHM for the broad neutral distribution 
relative to the one for the protons. HM07 further suggested that the FWHM versus Vs relation 
might be modified substantially if one takes the structure of the shock transition zone into 
account. If the velocity difference between the atoms and the protons were considerably 
less than the 3us/4 assumed in earlier models, then two consequences would result: 1. The 
amount of broad emission produced would be under-estimated; 2. For a given FWHM, 
the shock velocity inferred would always be less than the true value. However, as we saw 
in HSl the velocity of the proton fluid is decelerated to fH/4 almost immediately, and the 
subsequent evolution of v is weak. This validates the thin shock assumption made by CKR80 
and HM07 and highlights a puzzle — how does one account for the excessive amount of 
broad ("interior") emission observed in SNR 1987A (Heng et al. 2006; HM07)? More 
optical/ultraviolet spectroscopic studies of SNR 1987A are needed to shed light on the issue. 

Another relevant quantity to examine is the spatial shift between the centroid^ of ^f, 
and S,n- In Fig. [6l we compute and ^„ for /ion = 0.5 and = 0.1 cm~^, over the range 
1000 < f H < 10, 000 km s~^, as well as for (3 = 0.1 and 1. (Again, as a matter of illustration, 
we assume Case B conditions for the narrow Ha line.) Then, for 0.1 < /ion ^ 0.9, we 
determine the dimensionless spatial shift, 0shift (Fig. [7]), which is the spatial shift normalized 
by L. At any given shock velocity, Oshift decreases with increasing /ion- This is explained by 
the fact that at high neutral fractions (low /ion), there are initially only a small number of 
ions available; the system drifts along until there are enough ions to create the broad neutrals. 
Hence, the length scale for the creation of broad neutrals is relatively larger, corresponding 
to a greater shift in the centroid of ^b- For a fixed value of /ion, ©shift decreases as the 
shock velocity increases, for v^i > 1000 km s~^. This is because charge transfer reactions 
are favored at lower velocities and the larger number of broad neutrals created ensures that 
the centroid of is shifted farther downstream. This behavior is not true for shocks with 
fH ^ 1000 km s~^, as the charge transfer cross section becomes nearly constant with velocity. 

Observationally, it should be possible to measure ©shift with WFC3 onboard the Hubble 
Space Telescope if one isolates the narrow Ha component with a narrow band filter. Such 
measurements can be used to constrain no in some SNRs. 



^For an arbitrary distribution F{x), the centroid is defined as the point Xc such that J^" F{x) dx = 
C^F{x)dx. 
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6.2. FUTURE WORK 

Though the velocity of the broad neutrals and ions has no spatial structure, the densities 
of the atoms, broad neutrals and ions vary across the width of the shock transition zone. As 
mentioned in §11 knowledge of the density structures is relevant to modeling Lya resonant 
scattering in young, core-collapse SNRs. Photons are produced within the transition zone 
and resonantly scatter with a path length /mfp ~ 10^^ cm tioyr T^/Joo < ^zone, where tioyr is 
the time since the supernova explosion in units of 10 years and Tej = 100 K Tgj^ioo is the 
temperature of the freely streaming ejecta. There is evidence for Lya resonant scattering in 
SNR 1987A (Michael et al. 2003; Heng et al. 2006); while we now know how the density 
evolves spatially, it is beyond the scope of this paper to model the scattering process and 
observed spectra. It is, however, worthy to note that for Vg > 1000 km s~^ ( eioss ^ 1) and 
/ion ^ 0.9, the approximate, two-component solution for nn = rjnUa becomes a good one. 
The interested reader is referred to Zheng & Miralda-Escude (2002) and Tasitsiomi (2006) 
for the physics of Lya resonant scattering. 

We have constructed a three-component model where we employed the simplifying as- 
sumption vb = V. At fH ^ 3000 km s~^, impact excitation and ionization become competi- 
tive with charge transfer, and the resulting skewed Maxwellian of the broad neutrals peaks 
at Vb > V. The quantity /q°° dx is sensitive to changes in vb, which is relevant to the 
determination of Lf,n. Therefore, a vb ^ v treatment is necessary (M. van Adelsberg et al. 
2007, in preparation). Before (CKR80 and HM07), the ratio of broad to narrow Ha raiej^ 
was assumed to be equal to the ratio of rate coefficients, as the densities of the atoms, ions 
and broad neutrals were assumed to be constant. Physically, a full three-component model 
will answer the following question: for two Balmer-dominated SNRs with the same shock 
velocity, do we expect the ratio of broad to narrow Ha rates to be the same if one is highly 
ionized and the other is largely neutral? 

There are several important aspects of the hydrogen emission from non-radiative shocks 
that remain to be explored. If a significant fraction of the energy dissipated by a shock goes 
into cosmic rays, a precursor will heat and accelerate the gas before it reaches the shock 
transition zone, altering the density and velocity jumps and changing the post-shock ion 
temperatures (Blandford & Eichler 1987; Drury et al. 2001). Evidence for such a precur- 
sor is found in the anomalously large widths of the narrow Ha lines in Balmer-dominated 
shocks (Sollerman et al. 2003). Furthermore, when a pre-shock atom is ionized upstream, 
it effectively becomes a "pick-up ion" , analogous to those observed in inter-planetary space 



•^In this paper, we define the cmissivities as functions of x, i.e., (,{x). The rate is defined as the value of 
^(x) at a fixed value of x. 
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(Kallenbach et al. 2000). These pick-up ions are preferentially accelerated to become anoma- 
lous cosmic rays. In SNR shocks, they may form an isotropic, mono-energetic population 
that might perturb the Balmer line profiles. 
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A. APPENDIX: RATE COEFFICIENT 

The reaction rate coefficient between a beam of atoms and the Maxwellian population 
of an ion species s {p for protons and e for electrons) takes the form: 

Rs{v,T;vn,(r,ms) =J / /m T; m^) 6{v2-vn + v) a {\vi - V2\) \vi - V2\ d^vi d^V2 



= 27r /m,o /o~ exp 



2kT 



a (Av) Av Vr dvr dvzi (Al) 



where d^Vi — 27TVrdVrdVz and Av = y v"^ + (vz — \vh — "^1)^- The velocity difference between 
the atoms and the centroid of the ions is |fH — "^1, which is 3t'H/4 in the models of CKR80 
and HM07. The coefficient in front of the MaxweUian is /m,o = {ms/27rkTf^'^. 

The rate coefficient summed over species is 

R = Rp{v, T; vk, (Ji,p + (Jt,p, rup) + Re{v, (3T; vn, cri,e, rUe), (A2) 

since the atoms and protons participate in both ionization ( "7" ) and charge transfer ( "T" ) 
events, while the atoms and electrons only interact via the former process (if one neglects 
impact excitation). With the exception of that for charge transfer to an excited state (Barnett 
1990), all of the cross sections are taken from Janev & Smith (1993). 



B. APPENDIX: TRANSFER FUNCTION & DRIFT LENGTH 

The probability for destroying a broad neutral via impact ionization is described by 

^ = -I^P. (Bl) 
dx Vd 
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Upon its creation at xq, a broad neutral has a Pq = Poixo) chance of survivaL The probabihty 
of it drifting to a location x is given by integrating the previous equation from xq to x: 



P{x,xo) = Po{xo) 



Ci + C2 exp (C3X0) 
_ Ci + C2 exp (Cgx) 



(B2) 



where a = RiBVn/Rvd, Ci = bu/au, C2 = 1 — Ci and C3 = a^/L. Since the mass flux of 
broad neutrals has to be conserved, we require 



POO 

/ P{x, Xq) dx = 1. 
Jxq 



(B3) 



If we solve for Pq analytically, it contains the hypergeometric function of Gauss, 2F1] we 
choose instead to seek numerical solutions of Pq. 



If we express the transfer function in the form 

Pix,Xo) = Poixo) exp 

then the drift length is 

Idix, Xq) = {x — Xq) I « lu 



X — Xq 

ld{x,Xo) 



Ci + C2 exp {C^x) 



(B4) 



Ci + C2 exp (C3X0) 



(B5) 



We also note that for /ion ^ 0.9, Id^x^xo) ~ fd/n(xo)-RiB is an excellent approximation for 
the drift length. Examples of Vd are shown in Fig. [HI 
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Fig. 1. — Schematic diagram of the shock transition zone, in the case of a strong shock. The 
width of the transition zone is on the order of the mean free path of interactions (charge 
transfer and ionization). The velocity of the ions goes down to 1/4 of its pre-shock value 
almost immediately, according to the Rankine-Hugoniot jump condition. The ion density 
first jumps by a factor of 4 to conserve momentum, then eventually evolves to a value which 
depends on the pre-shock ion fraction, /ion- 
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Fig. 2. — Left: Rate coefficients for conversion of hydrogen atoms to protons (p) and electrons 
(e), used in the two- component model. The velocity difference between the atoms and ions is 
Svb./A. As an illustration, we display the individual rate coefficients for charge transfer {Ren, 
atoms and protons) and impact ionization {Rih', atoms, electrons and protons) for /3 = 1, but 
only the total rate coefficient for /3 = 0.1. Right: Rate coefficients for interactions between 
atoms (H), ions (i) and broad neutrals (B), used in the three-component model. 
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Fig. 3. — Spatial variation of rjn, rj, rj^ and e in the three-component model, for various values 
of the shock velocity, vs = fn/lOOO km s~^, and /ion — 0.5. Like in the two-component model, 
u ^ 1/4 throughout the shock transition zone. 
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Fig. 4. — Rate coefficients for the narrow (-RHQ,n) and broad {Rna,ho + -Rna.bJ Ha lines, 
assuming Case B and A conditions, respectively. 
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Fig. 5. — Emissivity of the broad and narrow Ha lines from a shock with parameters repre- 
sentative of SN 1006 (see text). 
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Fig. 6. — Emissivity of the broad and narrow Ha lines and their ratios, for various values of 
the shock velocity, vs — i'h/1000 km s~^, /ion = 0.5 and Ua — 0.1 cm~^. 
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Fig. 7. — Dimensionless spatial shift, ©shift, between the centroids of the broad and narrow 
Ha hne emissivity profiles. 



-25- 




Fig. 8. — Examples of the drift velocity, Vd, for (3 — 0.1 and 1. The lower and upper bounds 
are shown for 1/4 < Vd/vH < 1. 



